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Scalable arrangements of nitrogen vacancy cen- 
ters (NV) in diamond remain an open key chal- 
lenge on the way to efficient quantum informa- 
tion processing, quantum simulation and mag- 
netic sensing applications at the quantum limit. 
Although technologies based on implanting NV 
centers in bulk diamond [1] or hybrid device ap- 
proaches [2] have been developed, they are limited 
in the achievable spatial resolution [1] or by the 
intricate complexity, respectively. Here we pro- 
vide an alternative solution for creating a scalable 
system of individually addressable NV centers 
based on the self-assembling capabilities of bio- 
logical systems [3, 4]. By using surface function- 
alized nanodiamonds we propose a new avenue to 
bridging the bio-nano interface. Taking benefit of 
the outstanding nanometer resolution of the bio 
self-assembling techniques together with the con- 
trolled creation even of 3-D spatial structures [5] 
paves the way towards numerous multiqubit ap- 
plications. We provide a detailed theoretical anal- 
ysis on the feasibility of multiqubit quantum op- 
erations in one and two dimensional nanodiamond 
arrays, exploiting the significant dipolar coupling 
on the nanometer scale and address the problems 
of decoherence, imperfect couplings and the ran- 
domness of the relative orientations of the NV 
center symmetry axes. We show that our scheme 
allows for the high-fidelity creation of entangle- 
ment, cluster states and quantum simulation ap- 
plications. In addition we present first exper- 
imental demonstrations of interconnecting nan- 
odiamonds using biological protein complexes. 

Nitrogen-vacancy (NV) centers in diamond form 
promising candidates for quantum applications with the 
invaluable advantage of operating at ambient conditions. 
Ultralong electron spin coherence times in the millisecond 
range at room temperature [6] together with the ability 
of high fidelity polarization and optical detection of sin- 
gle electron and nuclear spins [7] allow for the coherent 
manipulation and read-out of the spin qubit states in this 
system. Interaction between NV centers can be achieved 



using the magnetic dipolar coupling [8, 9] between the 
electron spins of adjacent centers, yet this requires near- 
est neighbour distances of the order of 10 nm or below 
to enable significant coherent interactions. However the 
ability of current techniques (e.g. ion implantation) to 
position NV centers in very pure type Ha diamonds are 
limited to several tens of nanometers even for shallow im- 
plantation depths by fundamental straggling effects [10]. 

In contrast, the ability of biological systems for struc- 
tural self-assembly [3, 4] is a powerful tool allowing for 
the simple and parallel creation of large ordered arrays on 
nanometer scales thus outperforming the limited resolu- 
tion and structural complexity achievable of conventional 
lithography based on serial pattern creation. We pro- 
pose combining the spontaneous self-assembly of biologi- 
cal systems with the attachment of surface functionalized 
nanodiamonds. The biomolecules are used as a structural 
scaffold to enable the formation of NV center configura- 
tions with high spatial resolution. This can be performed 
by the combination of tiled motifs such as short DNA 
strands [11] or membrane forming complexes including 
SP1 [12], LH1 [13] and TF55/3 [14], that allow for the cre- 
ation of one and two dimensional arrays and have been 
successfully combined with nanoparticles. Going beyond 
two-dimensional periodic patterns the method of DNA- 
origami [5, 15], based on folding a large single stranded 
DNA molecule directed by staple strands, enables the 
creation of more complex highly controllable structures, 
ranging from aperiodic arrays to real three dimensional 
structures with nanometer precision. Elaborate knowl- 
edge in genetic engineering now allows control over the 
assembly, structure and topology of biomolecular com- 
plexes as well as the attachment of nanoparticles [16]. 
Each of these structures are suitable for hosting nan- 
odiamonds whose chemical attachment and site-specific 
binding can be directed by surface functionalization and 
labelling, and whose size can be controlled down to 4 nm 
in fabrication [17]. This allows for the high precision po- 
sitioning of NV centers interacting by dipolar couplings 
and paves the way for a highly controllable array of in- 
teracting quantum systems. 

Formation of SP1 nanodiamond complexes — 
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To demonstrate the feasibility of interconnecting nanodi- 
amonds with biological structures, we present the con- 
trollable and repeatable formation of small nanodiamond 
complexes using an SP1 (Stable Protein 1, see Supple- 
ment sect. 1) protein variant and first steps towards the 
formation of regular arrays of nanodiamonds on SP1 ar- 
rays. 

A crucial first step to enable both experiments is the 
genetic modification of SP1 to fuse 12 graphite-specific 
binding peptides to the SP1 N-terminus that permit site 
specific binding of SP1 to the carbon sp 2 -hybridization 
that forms on the nanodiamond surface [19]. For nanodi- 
amond samples (average size 30 nm) whose size exceed 
the diameter of SP1 (size 11 nm), small diamond clusters 
form when mixing solutions of nanodiamond and SP1 
performed under ambient conditions. The average num- 
ber of nanodiamonds in a complex can be controlled by 
the concentration ratio of SP1 to nanodiamond. Mixing 
a lmg/ml nanodiamond solution with a lmg/ml SP1 so- 
lution at a ratio of 1:1 mainly leads to the formation of 
30nm nanodiamond dimers and trimers (fig. 1 b-d) . As 
the nanodiamonds are larger than the SP1 the exact 
structure of the clusters is controlled by the nanodia- 
mond shape as several SP1 will bind the nanodiamonds 
across a surface. To verify that the creation of the clus- 
ters is due to the binding with SP1 rather than electro- 
static forces we have also followed the same procedure in 
the absence of SP1. In this case no nanodiamond aggre- 
gates are observed (fig. la). Site-specific binding due to 
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FIG. 1. SP1 nanodiamond clusters SEM images of small 
nanodiamonds complexes, using SP1 proteins to connect the 
nanodiamonds to each other, (a) no SP1 added gives single 
nanodiamonds. (b),(c),(d) mix of lmg/ml nanodiamond so- 
lution with lmg/ml SP1 solution at a ratio of 1:1 gives dimers 
and trimers. 

genetic engineering of the SP1 is essential for small nan- 
odiamonds (under 10 nm) to form regular structures on 
an SP1 array [20]. In a second experiment towards the 
creation of large ordered diamond structures, we achieve 
first examples of a dense nanodiamond monolayer (us- 
ing laser ablation made 5nm nanodiamonds, see fig. 2). 
Here a monolayer of genetically modified SP1 was formed 
by the Langmuir-Blodgett method [21] and subsequently 
combined with the nanoparticle solution. Note that the 
limited imaging resolution does not yet allow to verify the 
ordered structure of the nanodiamond array. Currently 



salts from the solution used in the Langmuir-Blodgett 
method crystallize on the substrate, leading to lattice de- 
fects in the nanoparticle array. These salt crystals may 
be removed by intensive washing, however with the risk 
of damaging the initial SP1 array arrangement (see Sup- 
plement, sect. 1). 
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FIG. 2. Towards ordered nanodiamond arrays (a) 

SEM image of a monolayer imaged after deposition of nanodi- 
amonds on an SP1 monolayer, (b) Schematic of an ordered 
hexagonal array of SPl-nanodiamond hybrids consisting of a 
nanodiamond attached to the SP1 inner cavity. 

Decoherence and dynamical decoupling - A 

key challenge for quantum computation and simulation 
with self- assembled arrays of nanodiamonds are the cur- 
rently achieved relatively short coherence times of a 
few lis [17] compared to the usual coupling strength of 
several kHz between NV-centers in different nanodia- 
monds [8, 9]. This is originated in interactions of the 
NV-center with the very proximate surface spins [22] and 
external charges, such that substantial improvements can 
be obtained by means of surface functionalization meth- 
ods [17, 23]. The remaining decoherence after optimiza- 
tion of material design can be further improved by ap- 
plying dynamical decoupling techniques, that allow to in- 
crease the coherence time by several orders of magnitude. 
These may be implemented either as pulsed schemes [24] 
or by continuous driving fields [9, 25, 26]. 
The pure dephasing noise can be modelled as a fluctuat- 
ing magnetic field energy shift b(t) with zero mean and 
(b(t) b(0)) = b 2 e~ l l T (see Supplement section 5). Adding 
a continuous driving interaction for decoupling, i.e. a 
resonant microwave coupling Vt on the relevant electron 
spin transition, the Hamiltonian of the effective two level 
system in a frame rotating with the laser frequency for a 
single NV center is given by H = (b(t)/2) a z -\- (Q/2) cr x , 
where a XjZ are the usual spin- 1/2 Pauli operators. In 
the relevant Markovian limit t ^> r and Qt > 1 the 
decoherence decay rate R(t) is determined by the noise 
spectrum S(u) evaluated at the decoupling frequency, i.e. 
R(t) oc S(Q) and the more general formalism will be dis- 
cussed in the Supplement (section 2) [27]. This allows to 
define an effective T2 time of the system resulting in 
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for the considered Lorentzian noise spectrum. Fig- 
ure 3 (b) compares this formula to numerical noise Simula- 
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tions showing that both the T2-scaling and the exponen- 
tial decay behaviour are in good agreement with numeri- 
cal calculations. Here the noise spectrum was calculated 
for a fluorine terminated surface of a spherical diamond 
with radius r = hum that results in fluorine nuclear- 
1/2 surface spins with a nearest neighbour distance of 
2.5 A, using the mean-field approach described in [28] (see 
Supplement, section 4). We thus expect a noise correla- 
tion time of r = 2.5/is and a mean square root ampli- 
tude b = 30.2 kHz (r ex l/(n 3 / 2 ) and b 2 oc n 2 r with n 
the surface spin density and r the nanodiamond radius), 
that could be further improved by spin bath polariza- 
tion schemes. For those parameters, a decoupling field 
of Q = 1.2 MHz leads to an effective coherence time 
of T2 c± 4 ms comparable to the electron spin Ti-time 
that imposes an upper limit on the dynamical decoupling. 
In experiments it is possible to achieve much higher de- 
coupling fields up to 300 MHz such that even a small 
number of possible electron spins can be decoupled (for 
a discussion of electron spin noise see Supplement, sec- 
tion 4). Combining self-assembled structures with litho- 
graphic techniques might help to integrate current struc- 
tures within the system [29] beneficial especially for high 
driving fields. Importantly, decoherence due to inten- 
sity fluctuations of the decoupling field can be suppressed 
by using a concatenated decoupling scheme as proposed 
in [26]. 
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FIG. 3. NV center, dynamical decoupling and dipolar 
coupling strength (a) Ground state electron spin triplet of 
the nitrogen vacancy center with zero field splitting 2.87 GHz 
and the degeneracy of the | =b 1) states lifted by a weak mag- 
netic field, (b) Effective T2 time obtained by dynamical de- 
coupling vs strength of the decoupling field: The continuous 
line corresponds to formula (1), circles correspond to the 1/e 
decay time and squares to an exponential fit both obtained by 
a numerical noise simulation. Noise parameters follow from 
the noise spectrum (inset) calculated using the method pro- 
vided in [28]: r 2.5 fis, b 30.2 kHz, T 2 (Q = 0) = 13.3 /xs. 
(c) Dipolar coupling parameter £ vs the external magnetic 
field. Red circles denote the average and blue lines the vari- 
ance arising from the random axis orientation. For large mag- 
netic field the quantization axis is given by the external field 
and £ ~ 1/4. The green dashed area indicates the range that 
can be corrected using compensation methods. All parame- 
ters are obtained by averaging over 10 6 random spin orienta- 
tions. Inset: Optimal magnetic field configuration for a 2D 
array also used for the main graph. 



Engineering of interactions, spin gates and 
quantum simulation — Dipolar interactions between 
electron spins of two adjacent NV-centers provide a pos- 
sibility for implementing gate operations [9] . Combined 
with the continuous driving of the decoupling field, the 
total effective Hamiltonian in a frame rotating with the 
microwave frequency can be written as (see Supplement, 
sect. 3) 
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Herein the last part accounts for the dipolar coupling 
with Jij = 2£ijHij where ^ depends on the external 
magnetic field strength] B | and its orientation with re- 
spect to both the NV symmetry axes and the vector 
fij connecting NV center i and j. In the limit of high 
magnetic fields ^ = 1/4(1 — 3 cos 2 ^), where Oij = 
Z(rij,B) and /i^- = fioj^h/ (Airr^) with the latter being 
Hij =52 kHz for an NV center distance nj — 10 nm. /xq 
and j e i denote the magnetic permeability and the elec- 
tron gyromagnetic ratio, respectively. It is important 
to note that the configuration of nanodiamonds leads 
to a random relative orientation of the symmetry axes 
in space that has two important consequences: In or- 
der to obtain a uniform dipolar interaction the quan- 
tization axis must be determined by the external mag- 
netic field as illustrated in fig. 3 (c) together with the 
optimal 2D-orientation, which requires that B > 0.5 T 
(7 e B /(2tt) > 14GHz), such that the magnetic energy 
shift j e B > D with D ~ 2.87 GHz the crystal field 
dependent zero-field splitting. Furthermore the orienta- 
tion dependent zero-field contribution leads to different 
transition frequencies of the individual NV centers with 
typical magnitudes differing by several hundred MHz, 
thus providing the possibility for individual microwave 
addressing for the typical values of Q obtained from fig. 3 
and at the same time maintaining a uniform dipolar cou- 
pling interaction. 

Combining the dipolar interaction with decoupling sup- 
presses the environmental coupling, i.e. decoherence, 
but also part of the gate interaction, a general prob- 
lem in the application of decoupling sequences. Trans- 
forming to an interaction picture with respect to the 
driving in the relevant limit Q ^> , the effective 
interaction up to non-nearest neighbours is given by 
Hi Mk = 1/2 V : ,, : ./„.^ u . S« i = s%s j _ +h.c. and 

^M 2 ~ s \ s \ + h ,c - wnerem (hj) sums over all neigh- 
bouring NV centers and Aii corresponds to the situation 
of a homogeneous decoupling field Qi = Q whereas A^2 
corresponds to the situation 1^ = — Qj for i G neighb(j). 
5 + and S- are the ladder operators in the a^-eigenbasis. 

The two qubit setting is shown in fig. 4 (a) (b), illus- 
trating the fidelity and time evolution to create a maxi- 
mally entangled state by a tt/2 and 57r/2 two qubit rota- 
tion in the Mi - coupling manifold, respectively, depicted 
in the left inset of figure 4 (a). Note that the required 
time for the latter case exceeds T2(£l = 0) by more than 
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a factor of three therefore impressively demonstrating the 
decoupling that allows for a 98% final state fidelity for a 
decoupling field of tt = 1.2 MHz. 
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FIG. 4. Two qubit gates (a) Fidelity vs decoupling field 
for a two qubit gate interaction and a tt/2 and 5tt/2 rotation 
in the Mi manifold leading to a maximally entangled state. 
Inset: Energy levels and dipolar coupling for the different 
manifolds, (b) State population vs time for zero decoupling 
(dashed) and Q — 1.2 MHz. The former case leads to a max- 
imally mixed state, (c) Fidelity vs systematic error e for a 
7r/2 rotation with (red) and without (blue) applying the error 
compensation sequence for a decoupling field Q — 0.5 MHz 
(continuous), Q = 0.8 MHz (dashed) and Q = 1.2 MHz 
(dashed-dotted). The noise parameters are given in fig. 3 and 
J = J12 — 26 kHz corresponding to the configuration in fig- 
ure 3(c) for ri2 = lOnra. 
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FIG. 5. Multi qubit gates (Noise parameters see fig. 3 
and Jij — J — 26 kHz corresponding to the configuration 
in figure 3 (c) for nj = 10 nm) (a) Fidelity for creating a 
one dimensional cluster state vs decoupling field strength and 
the number of qubits involved. Time addition is performed 
by two Suzuki- Trotter cycles. The inset shows the number 
scaling for a constant decoupling field of Q, = 0.5 MHz (blue 
continuous), Q = 0.8 MHz (red dashed) and Q = 1.2 MHz 
(green dashed-dotted). (b) Four qubit 2D cluster state for one 
(red) and two Suzuki- Trotter cycles (blue) eliminating next 
nearest neighbour interactions using the addition sequence 
sketched below: Herein blue lines denote a z a z interactions, 
red lines Mi interactions and dashed lines denote a negative 
sign of the corresponding interactions obtained by red circles 
representing embedding the interaction in between a U, 
pulse with U = exp(— i tt/2 a x ) on the red marked qubits. 
The basic interactions are obtained by adding Mi and M2 
as described in the text. 



In contrast to these simple two qubit examples, re- 
covering the full dipolar interaction form (<j z (g) <j 2 -type 
coupling, see eqn. (2)) is crucial for a wide range of appli- 
cations ranging from cluster state computation to quan- 
tum simulations. Once achieved this allows to create 
all other types of interactions by merely applying local 
unitary transformations. The <j*cr^-type interaction can- 
not be recovered by any local operation out of the re- 
duced manifold interactions Hi^ h ; however one can add 
up the contributions Hj ^i an d Hj m 2 m time by using 
the Trotter or Suzuki- Trotter formalism to implement the 
total Hamiltonian (see Supplement, sect. 6): 

H zz = H IMl + H IM2 = ^(Jij/2) *iai . (3) 

(id) 

For this scheme to work, it is crucial that the time ad- 
dition is based on the interaction frame Hamiltonians 
Hj^Mk instead of the one defined in (2), because only in 
that case the different timescales of the decoupling and 
coupling strength allow to create a high fidelity a\ a{ gate 
at the same time preserving the decoherence decoupling 
effect (see Supplement section 6.1 for more details). 

Compensation of systematic errors: Distance 
variations between adjacent vacancy centers and differ- 
ent orientations of the symmetry axes make it hard in 
practice to guarantee a uniform coupling. Therefore the 
coupling coefficient appearing in Hj^ k (for the two 
qubit situation) and H zz may be replaced by J(l + e^) 
with €ij describing the systematic error from the optimal 



case. However, extending the concepts provided in [30] 
allows to construct compensation sequences (see Supple- 
ment, sect. 7) provided that the error e^- < 1 and that 
non-nearest neighbour couplings can be efficiently sup- 
pressed, as can for example be achieved by approaches 
like the one discussed in fig. 5 (b). We analyzed the com- 
pensation method for the two qubit case in fig. 4 (c) and 
applied it to a four qubit cluster state in fig. S9 (Supple- 
ment). Due to the significantly increased time, that is 
eight and sixteen times the original gate operation, re- 
spectively, the region of benefit increases with the decou- 
pling field strength provided that it exceeds a threshold 
magnitude. As expected the two qubit gate compensa- 
tion is more efficient providing good results already for 
a 1.2 MHz decoupling field in contrast to the multiqubit 
counterpart that relies on a more general sequence less 
efficient in time. 

Cluster-state creation and quantum simulation: 
Once it is feasible to create a g\ (^-interaction, a cluster 
state can be realised, that subsequently allows to per- 
form quantum computations by purely local measure- 
ments. Provided that = J, cluster states can be cre- 
ated by a 7r/2-pulse of the o\ a 3 z - coupling, i.e. \cj) K=0 )c = 

S\ + + • • • +) with [31] S ~ exp (-i tt/4 £ ( . ?j) a\ trj) up 

to local operations that might be implemented by fast 
pulses on the electron spin manifold. Therefore the only 
complication consists in creating the Ising-type interac- 
tion by using the methods described in the previous para- 
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graph. As can be seen in figure 5 (a) and (b) the com- 
bination of decoupling and time addition is capable of 
creating one and two dimensional cluster states with fi- 
delities well above 90%. As a second application the sim- 
ulation of a Heisenberg-chain Hamiltonian is illustrated 
in the Supplement fig. S8. 

Summary — In summary we have demonstrated a 
new method to scale up quantum systems based on NV 
centers by exploiting the ability of biological systems 
for self-assembly with the precise positioning of nanodi- 
amonds in such structures. We experimentally realized 
the creation of small nanodiamond clusters as well as first 
steps towards large ordered arrays. Based on the achiev- 
able NV distances on the nanometer scale we theoret- 
ically proposed the realization of single and multiqubit 
gates as well as interesting applications such as the cre- 
ation of cluster states, addressing the typical problems 
as the limited coherence time and heterogeneous dipo- 
lar coupling strengths. Decoupling fields around 1MHz, 
well within reach in current experimental setups, to- 
gether with significant dipolar couplings of several tens 
of kHz, allow for the efficient decoupling of nanodiamond 
surface noise leading to coherence times comparable to 
the ultimate T\ limit. Thus gate fidelities well above 95 % 
can be expected even for multiple qubits and imperfect 
couplings. We believe that the combination of nanodia- 
monds with biological systems provides a promising ap- 
proach towards scalability overcoming the limitations of 
current attempts and offering a high level of control in 
both structure formation and qubit addressing. 
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1. SP1 PROTEIN COMPLEX AND EXPERIMENTAL DETAILS 

One remarkable aspect of biomolecules is their ability to recognize a wide range of substances with a high degree 
of specificity This feature of biomolecules has been extensively explored and a variety of artificial peptide aptamers, 
that specifically recognize various inorganic materials, have been created by selecting binders from random arrays of 
amino acids displayed on phages or bacteria (combinatorial biological method) [1—5]. 

SP1 is a thermally stable protein, originally isolated from poplar trees [6], which self-assembles to an 11 nm ring- 
shape dodecamer (12-mer). The protein is exceptionally stable under extreme conditions, being resistant to proteolysis, 
high temperatures, organic solvents and high levels of ionic detergent [7, 8]. SP1 multivalency allows the display of 
12 binding sites upon one protein complex, resulting in a stable and strong binding agent. 

Fusion of Carbon Nano Tube (CNT) specific binding peptides to SP1 N-terminus by genetic engineering resulted 
in an SP1 ring with 12 CNT binding sites, 6 on each side of the ring. This enabled the creation of SP1 variants which 
tightly bind to CNTs to form a stable SP1/CNT complex [9]. In this work we used the same SP1 variant to attach 
and order the nanodiamond structures. The carbon sp 2 -hybridization formed on the surface of the nanodiamonds 
was used to link the nanodiamonds. 

The formation of the small nanodiamond complexes was achieved by mixing 50 fiL of 1 mg/mL of the SP1 solution 
with 50 \iL of 1 mg/mL of the 30 nm nanodiamonds solution. To the mixture we added 900 jiL of distilled water. 
10 \iL of the final solution was applied on a silicon chip and scanned by scanning electron microscopy (Extra High 
Resolution Scanning Electron Microscopy Magellan™ 400L) as shown in figures 1 a-d in the main text. 




FIG. SI. (a) AFM scan and line profile (inset) of a dense monolayer of CNT binding SP1 with a scratched area, (b) Small 
ordered area of 5nm gold nanoparticles on an SP1 monolayer formed by the Langmuir-Blodgett method. The periodicity is 
observed to be 11 nm, equal to the diameter of the SP1. Alongside the nanoparticles we still have salt particles (the bigger 
particles), which damage the order and periodicity of the array. 

Formation of small ordered nanoparticles areas on an SP1 monolayer is performed using the Langmuir-Blodgett 
method [10]. In this method a trough is filled with a subphase solution that contains the SP1 proteins and by adding 
glucose to the subphase solution the SP1 floats to the solution-air interface. Then, by slowly reducing the interface 
surface area, the SP1 are forced to compress to a point where they become a monolayer. The SP1 monolayer is 
transferred to a silicon chip and washed with distilled water to remove excess salts from the subphase solution. In 
figure SI (a) an AFM scan with a scratched area is shown. By measuring the height difference between the silicon 
chip's surface (scratched area) to the rest, we verified the existence of a monolayer of height 2nm which corresponds 
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to the height of the SP1 protein on a Si surface. This shows that at this stage we have a dense monolayer of SP1 on 
the surface. In a second stage the silicon chip with the SP1 array is inserted into a beaker with nanoparticle solution 
on an orbital shaker for 10 minutes, and afterwards it is washed and dried. With gold particles (5 nm nanoparticles 
Sigma Aldrich, see fig. SI (b)) a small ordered area was already achieved. Excess salts remained, however, and formed 
salt crystals that lead to lattice defects in the monolayer. As for nanodiamonds (5 nm produced by laser ablation from 
Ray Techniques Ltd, see fig. 2 in the main text) the excess salts were sufficiently removed by cleaning but a periodic 
nanodiamond order was not observed yet. The samples were scanned by scanning electron microscopy (Extra High 
Resolution Scanning Electron Microscopy Magellan™ 400L) [11] as shown in figure SI and figure 2 in the main text. 



2. DYNAMICAL DECOUPLING AND THE FILTER- SPECTRUM OVERLAP APPROACH 



The effect of dynamical decoupling or more precise the decoherence decay rate of a decoupled system can be 
described in terms of a filter function overlap, representing the effect of the decoupling field, with the noise spec- 
trum [12-14] (see figure S2). This allows a very illustrative analysis and description of the working principles of 
dynamical decoupling methods, opens the possibility of measuring bath-coupling spectra by designing appropriate fil- 
ter functions and allows for the optimization of decoupling sequences. The only restriction arises from the assumption 
of the weak coupling limit, i.e. by treating the noise influence in a perturbative way. More practically this means that 
the coherence time must be large compared to the noise bath correlation time, a condition that -dependent on the 
type of noise- might not be fulfilled for small decoupling fields. However for the parameters considered in this work, 
in general T2 > r, and hence those limitations are not relevant. Moreover the expressions are exact in any limit for 
the free induction decay and pulsed decoupling schemes provided that the noise is Gaussian as can be easily checked 
comparing the final expression with the ones derived in a non-perturbative way in e.g. [15]. A summary of the, in 
many cases somewhat lengthy, formulas will be given in section 2.4. 



2.1. Decoherence decay rate 



Consider a two level system evolving under the Hamiltonian (in the rotating frame) 

H^K S f., + Kf„ (2.1) 

with S(t) a random, zero-mean fluctuating detuning describing the effect of pure dephasing and originating from 
environmental coupling and Q(t) the classical control decoupling field. We will assume that the system is initially 
(to = 0) prepared in the state 

m = +=(\e)+e i +\g)) (2.2) 

whose special cases = and = tt/2 correspond to the a x and a y (\+) y ) eigenstate, respectively. Then, 

after a time t the probability for still finding the system in that initial state is given by 

= I (l + a)s 2 0e-^W* + sin 2 0e- 1 /2(^W+fl7(*))*) , (2.3) 

describing purely the effect of decoherence and not taking the coherent evolution into account (see figure S3). The 
decay rates appearing in (2.3) can be expressed as 

1 1 f°° 

R k (t) = Yt ^ J 5 M dw ( 2 - 4 ) 

with S(lj) = ex.p(—iur)(6(t)S(t-\-r))dr the noise spectrum and F^(u) a decoupling field dependent filter function 
(see figure S2). Exact expressions can be found in section 2.4. 

R x (t) corresponds to the decay rate for an initial o x eigenstate with R x (t) ~ 1/2 S(Q) for t ^> r and a constant 
decoupling field, wherein S(Q) is the noise spectrum evaluated at the frequency of the decoupling field. The decay 
rate R y (t) = 1/2 (R x (t) + R 7 (t)) corresponds to the decay rate for an initial a^-eigenstate and is split into two parts. 
Splitting the decay rate R y (t) into the two contributions R x (t) and R 1 (t) allows for a more simple interpretation of 
the decay behaviour. In the limit of t ^> r and Qt ^> 1 the contribution R 1 (t) ~ (or more precise R 1 (t) t <C R x (t) t) 



S3 



and thus the decay rate of the second decay contribution corresponds to half the decay rate of the first contribution 
or as a specific example the decay rate of the a x is twice as large as the decay rate of the a y eigenstates. Therefore 
in that limit (t ^> r, Qt > 1) and for a constant decoupling field (2.3) takes the form: 

<Vvl/#0> = \ (l + cos^e-^^ + sinVe- 172 ^') with R x = ^S(Sl). (2.5) 

This has a clear interpretation in that the first part (the <j x eigenstate decay) is related to a population decay 
whereas the second part (the a y eigenstate decay) corresponds to a decay of the corresponding coherences, a process 
that happens in the Markovian limit t ^> r -also characterized by a time independent decay rate- with half of the 
population decay rate. In the non-Mar kovian limit the decay rate of the coherence contributions is increased by an 
additional factor R 7 (t). Note also that for the case of free induction decay, i.e. Q = 0, R 1 (t) = R x (t) and thus both 
decay contributions are equal as expected by the isotropy of the system. 
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FIG. S2. Decoupling effect illustrated in the filter - noise spectrum overlap approach for different times (a) t = 10 /is and (b) 
t = 100 /is. The red dashed curve corresponds to the noise spectrum whereas the blue and green curves represent the filter 
functions for a zero decoupling field and Q = OA MHz, respectively. For long times t ^> r the filter function becomes ^-peaked 
around uo — ±Q. Noise parameters: T2 — 13.3 /is, r = 2.5 /is, yj (5 2 ) = 30.2 kHz. 
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FIG. S3. Comparison of the decay rate according to (2.3) (red lines) to a numerical simulation of the noise process (blue) for 
an initial state (2.2) with (a) = 0, (b) = tt/2 and (c) <j> — 7r/4, Q = 0.5 MHz and modelling the noise as an Ornstein- 
Uhlenbeck process. Other lines correspond to decay curves for a zero decoupling field and different initial states as indicated 
in the figure. Noise parameters: T2 = 13.3 /is, r — 2.5 /is, \J (5 2 ) = 30.2 kHz. 



2.2. Decay rate derivation 



Following the description in [12, 13] we will derive the general decoherence decay formula (2.3). It is advantageous 
to evaluate the decoherence behaviour in the <j x eigenbasis = |d=) in which the effect of the decoupling field can 
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be interpreted as creating an energy gap suppressing flipping processes by the (off-resonant) detuning fluctuations. 
Using the Nakajima-Zwanzig projection operator approach [16] allows to obtain a master equation for the system 
interacting with a noise bath up to second order in the coupling constant for the evolution under (2.1) [12] 



dp 

dt 



= -Jp Jo dt ' ^ ~ ^ [S{t) ' S{t/) P{t)] + h '°' 



(2.6) 



with <f>(t - t') = (6(t)5(f)), S(t) = h/2a z and 

a z = exp (i J t Q(t')/2dt'a x ^ <x z exp ^— i J* Q(t')/2 dt' cr x ) the a z operator in the interaction picture with respect to 

the decoupling field. It turns out that coherences and populations are decoupled in the differential equation expressed 
in the er^-basis states leading to (pij = (i\p\j)), 



wherein (see definitions in section 2.4, U = exp(z J* fl(r) dr), 1Z denotes the real part) 

7l (t) = f dt' <t>(t -t')Tl [U(t) t/t (t')], R x (t) = \ f dt' 7l (f ') dt' 
Jo t Jo 



(2.7) 



(2.8) 



leading to the solutions 



p ++ {t)= L -[{2p\ + -l)e- R ^ 



-{t) = - [(2p°_ 



l)e 



-R x (t) t 



(2.9) 



On the other hand the differential equation for the coherences takes the form 



A (p-+ - p+- 

dt V P-+ + P+- 



1 / 71W+72W z(/ii(t)-/i 2 W) 

2 U(/ii(t) + /i 2 W) 71W-72W 



P-+ - p+- 
P-+ + p+- 



(2.10) 



with the additional definitions (herein 1Z denotes the real and X the imaginary part) 



!2{t)= f dt'<j>(t-t')1l[U(t)U(t')}, 
Jo 

m(t)= [ dt' <j>(t-t')i[u{t)u\t')} 

Jo 

p 2 (t)= I dt'<j>{t-t')l[U(t)U{t')} . 
Jo 



1 r* 

R 1 (t) = - / dt' 72 (t')dt' 
* Jo 



(2.11) 



As will be seen later, the combination p |_ — p^ is responsible for the decay of coherences in the density matrix 

description. Moreover it is straightforward to show that the off-diagonal elements (the //-terms) are related to a 
coherent evolution whereas the diagonal elements describe the decay terms of the corresponding quantities. Since we 
are not interested in the coherent evolution it is possible to set those off-diagonal contributions to zero resulting in a 
description of the envelope decay of the quantities. Therefore one ends up with 



^ (p-+ - P+-) = - \ [7i CO + 72 CO] (p-+ - P+- ) 



(2.12) 



(p_ + -p + -)(*)=exp(-± J\t'[ 7l (0+72(0)) (P-+ 

■p+-)(0). 



■P+-)(0) 



(2.13) 



Now consider the arbitrary phase state (2.2) that can be expressed in the |±) basis as (neglecting a global phase 
factor) 



|^)=cos(0/2)|+)-isin(0/2)|- 



(2.14) 
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leading to the initial density matrix 

Po = mm = l+X+l cos 2 (0/2) + l-X-l sin 2 (</>/2) 
-i(|->(+|-|+)<-|)cos(0/2) sin(0/2). 

Using eqn. (2.9) and (2.13), the density matrix after a time evolution of t follows to be 

Pit) = l+X+l \ [(2 cos 2 (</>/2) - l)e-«-W* + l 



(2.15) 



-)(- 



1 r 

2 



(2 sin 2 (</>/2) - l)e" 



-Rx(*)t 



(2.16) 



* (l-X+l - l+X-l) cos(0/2) sin(0/2) e -i/2(^(t)+K 7 (t))t _ 



Knowing the density matrix time evolution it is straightforward to calculate the decay rate of the initial state 
towards the completely mixed state pmixed — 1/2 (|+)(+| + |— )(— |) and the probability for finding the system in the 
initial state after a time t: 



tr(|^)<^| p(t)) = = \ (l + oos^e-^-W* + S in 2 0e- 1 /2(««(*)+^(*))*) 

what corresponds exactly to the result given in eqn. (2.3). 



(2.17) 



2.3. Population decay for = and — tt/2 and a constant decoupling field 

In this section the limiting cases of starting in the a x eigenstate \+) x (0 = 0) and the a y eigenstate \+) y (0 = 7r/2), 
respectively, shall be reviewed for the case of a constant decoupling field Q(t) = O = const., leading to an explanation 
of the different decay rates in both basis states. 

For = the initial state is given by (2.15) 

Po = l+X+l (2-18) 

and the time evolution follows from (2.16) 

^(*) = l+X+l I (l + e-^^ + l-X-l \ (l-e-^W*) (2.19) 

leading to a decay of the initial state (2.17) 

(+\p x (t)\+) = X - (l + e-^W) . (2.20) 

In the |d=) basis, that is optimal for the interpretation of the decoupling effect, this corresponds to a population decay 
with an effective rate determined by the flipping term S(t) (the pure dephasing in the original state basis) suppressed 
by the off resonance due to the decoupling field energy splitting Q. 
For = 7r/2 the initial state takes the form (2.15) 

pi = \ (l+X+l + l-X-l) - \ (l-X+l - l+X-l) • (2-21) 

Note that the first contribution already corresponds to the completely mixed state and does indeed not change in 
time such that the decay is determined by the decay of the a x basis coherences (2.16) 

P y (t) = \ (l+X+l + l-X-l) - \ (l-X+l - l+X-l) e"+ 2 (*+'+*+'))' (2.22) 

leading to the probability for finding the system in the initial state (2.17) 

tr {pi py{t)) = X - (l + e-+ 2 («*(')+«+')) *) . (2.23) 

Referring again to the decoupled decay picture it is clear that the decay of the a y eigenstates is governed by the decay 
of the a x coherences. 
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This interpretation of pure population (<fi = 0) and coherence ((f) = tt/2) decay allows especially for a simple 
explanation of the decay rate difference in the Markovian limit (t ^> r, Qt > 1): In that case the coherence decay 
induced by the population decay is just half that rate as also follows by noting that in that limit R 7 (t) ~ and 
R x (t) = 1/2 5(0) = R x = const. For a Markovian master equation description that property follows from the 
conservation of the density matrix trace as can be easily checked analysing the master equation (i.e. the trace of the 
right hand side should be zero): 



dp _ Rx 
dt ~ 4 



(2cr_pa + + 2cr + pa_ - {a + cr_,p} - {a_cr + ,p}) . 



(2.24) 



2.4. Decay rates, filter functions and definitions 



Decoherence rate 

Definition: 



R x (t) = - f dt' f dt" </>(t' - t") K (U(t') U\t")) 
t Jo Jo 

= - f dt' f dt"^{t f -t f, y 

t Jo Jo 

ds ^ drft(r)^ cos ^ drft(r)^ + sin ^ drft(r)j sin ^ drft(r) 



Ry(t) = - f dt' f dt" W - t") K (U(t') U(t")) 
t Jo Jo 

= - f df f dt"cf>(t f -t"y 

t Jo Jo 

)s I dr ft (r) j cos (j dr ft(r) j - sin M dr ft (r) j sin N dr ft (r) 



Filter functions 



with 



and 



/CO 
e- iwT ((5(t)5(t + T)}dT 
-CO 

jf e- iwt 'cos^ drfi(r)j dt' +^ e" ia "' sin ^jf drft(r) 



dt' 



J Q e ~ iUt ' C0S (^J Q dTfi(r)j dt' -^e- iw *'sin^ drft(r)j dt' 



Normalized filter functions F^ norm (cj) dto = 1 



F ,„ H= 1 FfH 

Z 7T t 



F^ norm (uj) = F?(u>) with A/" 7 = 2^ J* dt' 2 cos 2 N fi(r)d 



v 7 

Specific filter functions: 
• Free induction decay (fi = 0) 



I 7T t ZTTt 



lim F™ orm (w) = (J(w), lim ^(i) = lim i? 7 (i) = - 5(0) 



Continuous control field = const.) 



F t x (u)= 1 -t 2 



smc — - — £ + sine 



7 2 cos(Qt) (cos(Qt) — cos(cj £)) 



F , OTH= i ^ H 

Z 7T t 



Ft ~ nsm(2nt) FAuj) 



Markovian limit: t^> t, Sit > 1 



lim F^'" (to) = I [5{u -Sl) + 5{u + SI)] 



lim R x (t) = ls(Sl) 

£— )-00 Z 



lim RJt) = 
Spectrum for the Ornstein Uhlenbeck process 

_ 1*1 



(5(t)5(0)) = b 2 e 



5w = — tt^ 

1 + UJ Z T Z 
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3. HAMILTONIAN, PSEUDOSPIN 1/2 DESCRIPTION AND QUANTIZATION AXIS ADJUSTMENT 

USING AN EXTERNAL MAGNETIC FIELD 



Let us first consider a system of NV centers, namely the electron spin-1 ground state triplet manifold ( 3 A) coupled 
by a dipolar interaction: 



H — Ho + Hdip 

with the zero-field and external magnetic field contribution [17] 



(3.1) 



(3.2) 



Herein D denotes the orientation dependent zero-field splitting tensor, Si the spin-1 operators, B the external magnetic 
field and j e i the gyromagnetic ratio of the NV center electron spin. In the principal axis system defined by the NV 
symmetry axis, the zero-field splitting tensor D is diagonal and takes the form: 



Dj — diae 



1 



-D 



--D-E,-D 
3 ' 3 



(3.3) 



with D = 2.87 GHz and E introducing an additional coupling capable of lifting the | ± 1) state degeneracy. In 
nanodiamonds the strain dependent E contribution can be nonzero in contrast to the bulk diamond situation [18]. In 
the principal axis frame (denoted by the primed quantities) Hq can be rewritten as 



ff = 5> 



S' 



1 -2 

3 S ' 



e [s:?-sZ]+ lel B'-s> 



(3.4) 



providing a good choice for either small magnetic fields or in cases where the magnetic field is parallel to the NV 
center symmetry axis. An arbitrary choice of the reference frame (and therefore the quantization in that case), e.g. 
choosing a constant laboratory frame for multiple qubits with different symmetry axes, requires to rotate the tensor 
correspondingly such that it takes a different form for each of the NV's (unless they do have the same orientation). 
The dipolar coupling term has the form [19] 



H, 



dip 



i>j 



4-7T rf 



Si Sj 3 ^Si • G-ij^j ^Sj ' G-ij^j 



(3.5) 



with /i the magnetic permeability and and e^- the distance and unit direction vector between NV centers i and 
j, respectively. Assuming an equal quantization axis and imposing the secular approximation allows to approximate 
eqn. (3.5) by 



H d 



i>j 

= £ 



- ^ 7 ¥) (l-3cos 2 %) 
2 \ 4tt rf I v ' ' 



3 S^ Sj Si Sj 



i>j 



1 f MO 7ef ft 

47r rf- 



(3.6) 



(1-3 cos 2 %) [2 SI S z j - (Sf S x 3 + S? S] )] 



with Oij the angle between e^- and the quantization axis that is given by the external magnetic field direction for 
B ^> D. Note that the dipolar coupling in the present situation is rather small, e.g. (/xq J^i ^/(^ 7rr ?j)) = 27r • 52 kHz 



for 



10 nm and therefore small inhomogeneous broadening effects (e.g. different strain contributions) as well as 



different orientations of the NV center symmetry axis will essentially reduce (3.6) to 



3 cos 2 i7 -) S?S* 




(3.7) 



Creating an arrangement of nanodiamonds has the disadvantage that the symmetry axis of the NV center is not 
controllable, i.e. there exists no common quantization axis in a laboratory frame. Therefore, without or with a small 
magnetic field j e i B <C 1 the individual spins have an equal probability of pointing in any arbitrary spatial direction 
and therefore (3.5), evaluated in the principal axis symmetry frames of the centers involved, leads to a coupling that 
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depends as well on the relative orientation of the NV symmetry axes as on the one to the vector connecting the 
centers, therefore leading to a vast distribution of coupling frequencies inside a mutli-qubit array. However for a 
strong magnetic field (7^ B > D) the quantization axis is determined by the direction of the external magnetic field. 
Considering the limiting case 7^ B ^> D (B > IT) allows to rewrite the Hamiltonian as (neglecting off-resonant 
coupling terms of the zero-field splitting contribution) 

H ~Y. \ ( S -~l^) & (1 + 3 cos^tfi)) +3£ (l-coa(2ti i ))]+>y el BS z + H dip (3.8) 

with $i the angle between the magnetic field direction z and the symmetry axis of NV center i and Hdi P given 
by (3.5)-(3.7) with Oij the angle of the vector connecting i and j to the z-axis, i.e. the external magnetic field. Two 
important prerequisites are achieved that way: First the quantization axis is now completely determined by the 
external magnetic field and not by the symmetry axis any more such that the dipolar coupling will be homogeneous 
and independent of the individual orientations (for equal %). Second it will be very unlikely that two NV centers have 
the same energy transition frequency, i.e. there exists are strong orientation dependent (i?^) distribution of transition 
frequencies differing in the GHz range that allows for individual addressing and moreover in general, recalling that 
the mutual coupling is of the order of several tens of kHz, reduces the dipolar coupling to the form given in (3.7). 
In the more general case of intermediate magnetic fields the effect on the dipolar couplings is more complicated and 
cannot be simply described in the S x1 S yi S z basis manifold. We analyzed this case numerically below for the relevant 
quasi-particle two-level system. Effectively, recalling that there will be a continuous microwave driving to reduce 
decoherence, only a reduction to two states of the Spin-1 system is relevant, i.e. S| = | + + 1/2 (crj — tj) or 

Sj = I — 1| + 1/2 (crj + tj) for a zero magnetic field or a magnetic field parallel to the NV center symmetry axis, 
and therefore the two level reduction allows to rewrite (3.7) as 

H Ip S - E ( £ (1 " 3 cos2 %) \ K °*i T < T o]) . (3.9) 

Note that the single flipping terms can be incorporated in the energy Hamiltonian Hq, will be suppressed anyway when 
adding a sufficiently strong continuous decoupling driving field or can also be removed exactly by an echo sequence 
as discussed in section 6. 

Adding noise and an additional continuous driving for decoupling, the total Hamiltonian in the two level basis can 
now be written as 

rrTLS rrTLS 1 rrTLS 1 rrTLS . rrTLS /o i n\ 

H tot - H + H drive + H dip (d.lUJ 

with Hq LS the energy of the two selected dressed states of (3.2), i.e. obtained by diagonalizing (3.2) 

i 

and ujf = (D + j e i B)/2 for a zero magnetic field or a magnetic field aligned with the NV symmetry axis (and E=0) 
or uf = (1/8) [D (1 + 3 cos(2#i)) + 3E (1 - cos(2^))] + lei B/2 in case of a strong magnetic field j ei B > D. 
H^ LS describes the decoherence influence modelled as a pure dephasing noise 

H TLS = J2 b M at (3.12) 

i 

with bi(t) a fluctuating frequency as will be discussed insection5. 

H drive describes the continuous resonant microwave driving for decoupling on the dressed state two level system (3.11) 
and is given by 

i 

with Qi the Rabi frequency of the driving. 

The dipolar interaction in the secular approximation can be, noting again that only af cr| terms will in general survive, 
written as 



i>j 



13 



SIO 



where in the limiting cases of a strong magnetic field, or more general for equal quantization axes of i and j the 
parameter ^ is given by &j = 1/4(1 — 3cos 2 ^j). For general magnetic fields the parameter ranges are plotted in 
figure S4 and figure 3 (c) in the main text, showing that around B > 0.5 T are required for providing an almost uniform 
coupling; noting that the variance is even lower there is a high probability that a sufficiently uniform dipolar coupling 
is already obtained for B > 0.2 T. 

Finally, transforming to a rotating frame with respect to the laser frequency and neglecting counter- rotating terms, 
(3.10) then reads 
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(3.15) 
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In the attempt to achieve a uniform dipolar coupling by means of an external magnetic field the orientation of the 
magnetic field is crucial as can be seen from (3.7) and (3.14). Claiming that the interaction strength should be equal 
in all spatial directions the maximal dipolar coupling is achieved in a linear chain and in a two dimensional array 
when the magnetic field is parallel to the chain = —1/2) and orthogonal to the plane (&j = —1/4), respectively. 
In a three dimensional array this cannot be achieved, because for the only choice that provides equal interactions in 
all directions, the space diagonal, cos0^- = l/\/3 and therefore the interaction strength equals zero (performing an 
addition in time with non-equal couplings in all spatial directions can however circumvent that problem). As a last 
remark it should be noted that the combination of a strong magnetic field and the condition cos0^ = 1/y/S could 
be used to decouple the system from the dipolar interaction, e.g. after a cluster state is achieved and further local 
operations and measurements are desired subsequently. 
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FIG. S4. (a) Dressed state energy gap uj s in units of the zero field splitting D — 2.8 GHz for different magnetic field strengths. 
The two level system was chosen as the transition of neighbouring dressed state energy levels with the largest energy gap. Blue 
lines indicate the absolute deviations depending on the symmetry axis orientation, (b) Dipolar coupling parameter £ vs the 
external magnetic field parallel to the axis connecting the NV centers. The red circles denote the mean values and the blue 
lines indicate the absolute deviation (left) and variance (right) due to the random orientation of the NV symmetry axis. All 
values are obtained by averaging over 5 • 10 5 random axis orientations. 



4. NOISE SPECTRUM FOR A FLUORINE-TERMINATED NANODIAMOND SURFACE 

Terminating the nanodiamond surface by fluorine, oxygen or hydrogen / hydroxyl groups replaces the electron spins 
associated with the sp2 hybridized orbitals by much weaker nuclear spins [20] and therefore reduces the dipole-dipole 
interaction strength between surface spins by a factor of 10 -6 (as given by the square ratio of the corresponding 
gyromagnetic ratios). Additionally the coupling to the central spin, the electron spin of the NV center, reduces by a 
factor of 10 -3 . Recalling that the pure dephasing noise responsible for the decoherence mechanism can be described 
by a fluctuating magnetic field and that the strength and fluctuation rate depends on the (hyperfine) coupling to the 
central spin and the surface spin flip-flop rate, respectively, a significant improvement can be obtained by terminating 
the surface purely by nuclear spins. As an illustrative example we will concentrate on the fluorine terminated surface 
(nuclear spin 1/2), having the additional advantage of not affecting the charge state of the NV center and forming a 
lattice with nearest neighbour distance 2.5 A [21]. 

Such a coupled system can be described by 

H = H» V + H* f + H*i t + H hf (4.1) 
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with Hq V the NV center energy Hamiltonian as given in (3.2), Hq the magnetic field splitting of the surface nuclear 
spins, H*J[ t the dipole dipole coupling between surface nuclear spins and H^f the hyperfine coupling of the surface 
spins to the central spin (the NV center electron spin): 

Hh f ~ 2^ [j^ ~^f- J - 3 cos °i) S a z • 

Herein S denotes the spin-1 operator of the vacancy center, Gk the spin 1/2 operators of the surface nuclear spins, 
Tij the distance between spins i and j and T{ the one between surface spin i and the central spin and Oij and Oi the 
angles between the vector connecting the two coupled spins involved and the external magnetic field. As discussed 
in section 3 we will assume that the quantization axis of the NV center is essentially determined by the external 
magnetic field and not by the symmetry axis of the vacancy center. Alternatively one might assume that the external 
magnetic field is parallel to the NV symmetry axis. In the two-level approximation, used for the driven system in 
section 3, one can substitute S z —> 1/2 (s z ± 1) with s z the spin 1/2 operator of the quasi-spin. Herein we just 
retained secular contributions of the dipole dipole coupling and neglected flipping terms in the hyperfine interaction 
(that would be related to Tl) due to the large energy mismatch between surface spins and the NV electron spins. 
However, nevertheless those flip flop processes are very unlikely, second order processes might lead to a significant 
enhancement of the surface spin flips described by H-l t of the order of A? / A with A being the typical energy splitting 
of the NV levels [22-24]. However we do not expect those terms to be significant in the present proposals that requires 
magnetic fields B > 0.5 T what leads to mediated couplings of the order of A^/Q ~ 1 Hz and has to be compared 
to the direct nuclear nuclear dipole coupling, that, in the dense surface arrangement (2. 5 A) leads to next neighbour 
flipping rates of 6.8 kHz and is comparable to the mediated one for a nuclear spin distance of 5nm. Thus, in the 
dense bath at large magnetic fields considered here, the decoherence effect arises mainly from the flip-flop interaction 
of close nuclear spins which is to a good approximation not affected by second order hyperfine processes whereas the 
decoherence influence of the weakly interacting remote spins, affected by the mediated contribution, can be neglected 
compared to the first part (despite the fact that it leads to larger differences in the field). 



mutual electron spin coupling: 

CeLel = [ , 

V 4tt J 


electron spin - fluorine 
nuclear spin coupling: 
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mutual nuclear 
spin fluorine coupling: 


52 MHz (nm) 3 


74AkHz (nm) 3 


106.3 Hz (nm) 3 



TABLE SI. Magnitude for different coupling constants. Note that the coupling to the NV center corresponds to the electron 
spin coupling. 

Calculating the noise spectrum and the corresponding noise parameters is in general intractable for a large number 
of surface spins due to exponentially increasing computational resources. As early as in 1962 Klauder and Anderson 
showed by very general arguments that a Lorentzian noise distribution has to be expected in case of dipolar couplings to 
a spin bath [25]. Since then various approaches, mean-field and exact approaches in certain limits, have been invented 
to infer the noise properties and calculating the decoherence decay rate, ranging from the simple Ornstein-Uhlenbeck 
model [26] to cluster expansion methods [22, 24, 27, 28]. In here we used the mean-field method described in [15, 29] 
what corresponds to the lowest order of the cluster expansion methods, the so called pair-correlation approximation, 
corrected by a mean-field broadening using the method of moments [30] . In the pair-correlation approximation one 
assumes that each of the flipping processes in the dipolar coupling (4.2) is independent of all other processes, i.e. the 
Hilbert-space is substituted by a pair Hilbert space (ij) wherein (Hij follows directly from eqn. (4.2)) 

Htlt = Hij = %) with \ H i3 > H ki] = Vij, kl . (4.3) 

i>j (ij) 

This leads to a discrete spectrum of £-peaks at the different pair-induced transition frequencies and can be justified as 
long as correlations between those pairs can be neglected, that is as long as the evolution time from an initial thermal 
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state is short enough such that 1 — ex.-p(—qN 2 lip /N) is small[24] (with Nfu p the number of flipped bath spins during 
the considered time, q the number of nearest neighbours and N the total number of bath spins). Higher orders would 
lead to additional frequency peaks as well as couplings of the already existing peaks, i.e. a finite lifetime broadening 
what is taken into account by using a mean-field type approach based on the theory of moments. From eqn. (4.2) and 
(5.1) the operator for the effective field on the NV center follows to be 



E 



MO lei 7n h 
47T rf 



(1 -3 cos 2 Qi)G\ 



(4.4) 



and the noise spectrum can be calculated by (see section 2.1) 

r-OO 

S(oj) = dt (b(t) 6(0)) 



(4.5) 



where in good approximation the initial bath states can be assumed to be uncorrelated with the NV center and at 
room temperature given by 1/2^1 [26]. Following the calculation presented in [15] this leads to (neglecting static 
contributions that were also neglected in the decoherence discussion in section 2.1). 
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and cfij the mean field broadening used in replacing the original delta-peaks by Gaussian peaks leading to the same 
second moment as the one obtained by the moment theory 
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(4.8) 



The spectrum (4.6) was numerically evaluated by randomly distributing nuclear surface spins on a sphere with equal 
inter-spin distance and placing the NV center in the center of the sphere, i.e. for a sphere radius r = 5 nm a number of 
5026 (what corresponds to a distance of 2.5 A) nuclear spins was distributed on the surface (see figure S5). Assuming 
an equal quantization direction for as well the NV center and the nuclear spins (corresponding to a strong magnetic 
field or to a magnetic field parallel to the NV symmetry axis) the noise spectrum was calculated and subsequently 
fitted to a Lorentzian in order to obtain the noise amplitude b and correlation time r defined in (5.3). Depending on 
the number of surface spins an average over different random positions is performed (what however does not lead to 
a significant change for several thousand spins considered here). 

We performed the same calculation as well for different numbers of electron surface spins (see figure S6) on diamond 
with radius 5 nm. In that case the decoherence parameter are very poor, e.g. T 2 = 0.1/is, r = 21.7 ns and b = 3.6 MHz 
for 30 surface spins and getting worse with further increasing the number of spins. For such a situation decoupling 
fields in the Q ~ lr ~ GHz range are required, that, despite still allowing simple two qubit gates (provided the 
ft stability can be achieved), are intractable for multi-qubit applications that require a certain level of individual 
addressing. Of course it is questionable to what extend the mean-field spectrum calculation, initially based on a pair- 
approximation, retains its validity in the fast flipping regime of electron spins. To gain some insight in the electron 
spin case we performed some exact numerical T2 calculations for small electron spin numbers (note that for such small 
numbers the results obtained by 4.6 are not justified because the spectrum obtained that way differs significantly from 
a Lorentzian and is very sensitive to the position of the surface spins): For 10 surface spins T2 is already below 0.5 [is 
such that the correlations times obtained by the spectrum approach, even they might be not exact, seem realistic at 
least in the order of magnitude. Therefore it is crucial to avoid the presence of surface electron spins for performing 
reliable quantum gates. 
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FIG. S5. Surface spin de-phasing influence: The dephasing due to the nuclear surface spins (blue) is modelled by randomly 
placing nuclear spins on the nanodiamond surface with the diamond assumed to be a sphere of radius r = 5nm. The NV 
center is placed in the center of that sphere with the symmetry axis aligned with the external magnetic field that determines 
the quantization axis of the surface spins. 
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FIG. S6. Electron spin decoherence: (a) T2 (blue) and T2 SE (Hahn-echo, red) time for different numbers of electron surface 
spins calculated using the mean- field broadening method described in the text. A minimal number of ~ 30 spins is required for 
a Lorentzian-type noise spectrum and the possibility to consider mean values out of the random distribution. Inset: T2-time 
out of an exact numerical simulation for a small number of electron spins emphasizing that T2 times far below lfis are expected 
even for a small number of 10 surface spins, (b) Noise correlation time vs number of surface electron spins. Error bars indicate 
deviations from the mean value due to the random choice of the electron spin positions. Inset: Noise amplitude vs the number 
of electron surface spins, (c) Estimated scaling according to (2.5) of the effective T2 time with the decoupling field f2 for 30 
electron surface spins. The noise parameters are obtained from the mean-field approach and given by T2 = 0.1 {is, r = 21.7 ns, 
6 = 3.6 MHz. For those parameters decoupling fields in the GHz range are required which are challenging concerning intensity 
fluctuation, individual addressing and the validity of the rotating wave approximation. 



5. DEPHASING NOISE SIMULATION AND THE ORNSTEIN- UHLENBECK PROCESS 



Dephasing of the nitrogen vacancy center spin states is mainly caused by long range dipolar interactions with a 
spin-bath, consisting e.g. of substitutional nitrogen atoms (PI centers) and 13 C nuclear spins [31, 32] or in the case 
of nanodiamonds predominately of unpaired surface spins [18, 20, 33]. Herein a significant energy mismatch of the 
transition energies prevents flip-flop processes between the vacancy center and the bath spins limiting the influence to 
a pure dephasing process. However intra-bath flipping processes are not suppressed and occur on the noise correlation 
timescale r such that the NV center is influenced by a random configuration of the bath environment what can be 
modelled in the mean field approximation as a random magnetic field leading to the frequency shift b(t) and therefore 



H = b(t)S z . (5.1) 

Now assuming that the back action of the nitrogen spin on the large bath is small and that there is no coherence 
between different bath spin evolutions, one can model b(t) as a random Gaussian, Markovian and stationary process 
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what is also known as an Ornstein-Uhlenbeck process [32] and can be specified by the two time steady state correlation 



(b(t)b(0))=b 2 e X p(-\t\/T) 



(5.2) 



with b a measure of the coupling strength and r the noise correlation time, or alternatively by the expected Lorentzian 
noise spectrum [34] 
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For such a process the evolution of the random field is governed by the Langevin equation [35] 

db(t) 1 



dt 



b(t) + y/cT{t) 



(5.3) 



(5.4) 



with T(t) a Gaussian zero-mean noise [36] ((T(t)} = 0, (r(t)T(t / )) = 5(t — t')). The parameter c in (5.4) is related to the 
definitions in (5.2) by b 2 = cr/2 what follows by solving (b(t) 6(0)) by means of (5.4). Exact updating formulas can be 
obtained for the differential equation (5.4) and also for the time integral of b(t) given by y(t + dt) = y(t) + b(t) dt [35]: 



b(t + At) = b(t) ii + £ x m 
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(5.5) 
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with 



H — exp(— At/r) 

e x = (cr/2)(i-v 2 ) 

e y =cr 3 [AVr-2(l- M ) + l/2(l-^ 2 )] 
K =( C r 2 /2)(l-M) 2 



(5.6) 



and ni, n2 statistically independent unit normal random numbers. 

Simulation of the noise process: Simulating the evolution governed by a Hamiltonian of the form (3.15) can 
be easily performed by numerically solving the Schrodinger equation and using the exact updating formulas (5.5). In 
here we used a Suzuki- Trotter decomposition splitting the Hamiltonian in the noise H n and remaining part H ri such 
that the time evolution is given by (h = 1) 
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(5.7) 



and the noise part can easily be evaluated using (5.5). 



6. CREATING AN EFFECTIVE a z a z TYPE COUPLING 

A full Ising interaction, that is a a z <j z type coupling in the presence of decoherence is a crucial task towards the way 
to a universal set of gates, that, once achieved, easily allows to create all other types of interactions by just applying 
local unitary transformations, straightforward to implement by local microwave pulses on the electron spin transition, 
e.g. : af a* = U x of oj or of a v - = U y of oj with U x = exp (i tt/4 (o* + a J y )) and U y = exp (i tt/4 (a l x + a J x )) . 

Noting that (see also definitions in the main text) 



H zz = H IMl + H IM 2 = ( J ^j/ 2 ) , 

i>j ,iGneighb(j) 



(6.1) 



allows to use the time-addition methods to implement the H zz time evolution by either using the Suzuki- Trotter 
formalism 



exp(-iH zz t) = [R Ml (t/(2n)) R M , (t/n) R Ml (t/(2n))] n + 0{{Jijt/nf) 



(6.2) 
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or Trotter formalism 



exp(-iH zz t) = [R Ml (t/n) R M2 (t/n)] n + 0((Jijt/n) 



(6.3) 



with RM k (t) = exp(— i Hi y M k The switching frequency between the two types of interaction, Mi and depends 
crucially on the Hamiltonian timescales and to avoid a destruction of the decoupling effect, it is important to work in 
the interaction picture frame of the Hamiltonian Hj^M k , i- e - 



dip 



- 9 E J a s Mk 
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wherein Hj^ k 



Y,i^i/ 2(J i represents the decoupling field, H dip = 1/2 ^2 i>jieneighh ^) Jij <* % z °i the dipole-dipole 



n.c, d M2 



h.c. with s± the corresponding ladder operators in the a x 



coupling term and S 1 ^ = s l + s 
eigenbasis. 

Recalling the interaction picture definition, i.e. the connection of the interaction picture state \%jj(t))i with the 
original one \^(t)): 
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(6.5) 



it is obvious that an interaction picture time evolution can be created by the following pulse sequence: 



R Mk (t) ■= e~ iH ^ 1 = e +! ^ f e- ,fl ^ * = E Mk (t) U Mk (t) 



(6.6) 



with UM k (t) — exp(— iHM k t) an d EM k (t) = exp(z Hj^ t). There are essentially three options to create (6.6): 
Directly implementing the pulse sequence, adjusting the total time and making use of the different timescales in 
the Hamiltonian parts or by using an echo sequence. While the first method is obvious we will discuss the last two 
approaches in the following. 

Note also that a larger Trotter slicing parameter n does not necessarily increase the fidelity; the time interval must 
still be large enough such that Vt (t/n) > 2tt and therefore the off-resonant contributions average to zero and one really 
ends up with a proper manifold interaction in the interaction picture frame. For the typical parameters considered in 
this paper Jij ~ 26 kHz and Q ~ 1 MHz optimal values for t = 7r/(2 Jij) are given by n ~ 2, 3. 

Time adjustment: Noting the two different timescales in the Hamiltonian Hj^ k ^ i-e. the one of the dipolar 
coupling in the kHz range and the decoupling field strength in the MHz range it is possible to create the effective 
interaction (6.6) by merely adjusting the total time. This can be seen by noting that the switching process in the time 

addition of Aii and M 2 happens on the timescale of ~ j\- 1 \ On this timescale, what forms the typical timescale 
of t appearing in (6.6), the Hamiltonian contribution Hj^ k can be viewed as performing multiple 2tt pulses that are, 
up to a global phase, not relevant for the interaction. Therefore the pulse EM k (t) can be implemented on a much 
shorter time t' determined by ft -1 on which to a very good approximation EM k (t') ~ /7^ fc (t / ). t' can be determined 
by the following consideration (note that global phase contributions are neglected) 
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(6.7) 



with t' = (2tt- mod [fit, 2 tt]) /ft. 

Therefore, if one wants e.g. to implement the Suzuki- Trotter sequence 
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one can replace EM k (t/(2n)) ~ ^M fc (^n) w ^ n = (2?r — mod [Ot/(2n), 27r])/fi leading to 
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or equivalently the total time £ by 
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(6.9) 



(6.10) 



(6.11) 



resulting in a simple evolution under the original Hamiltonian HM k with an adjusted total gate time. 

Echo pulses: Another possibility of implementing (6.6) can be performed eliminating the additional Hj^ contri- 
bution by echo pulses. Out of (6.6) it follows immediately that 



e -iH A4k t = e -iH Mk t e -iH IjMk t^ 

The goal will be to cancel the additional Hj^ k contribution by means of a unitary (echo) 7r-pulse SV, i.e. 

S^e'* Hm * t/2 S 1 " e~ l HM k l l 2 = e~ { [Sn H ° M * S ^ ] t/2 e~ l ^ H ^ k si) t/2 e ~% H° Mk t/2 e _« Hl , A4k t/2 
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(6.13) 



thus leading to the conditions 



(i) S w Hi tMk Si = H It M k 

(ii) S^H° Mk St = -H° Mk 
(Hi) [H IMk ,H° Mk }=0. 

It is straightforward to verify that those conditions are fulfilled by performing a 7r-pulse in a z or a y , i.e. 



(6.14) 



S n = exp 



or S n = exp 



(6.15) 



Hence the interaction (6.6) can be implemented using the following sequence 

R Mk (t) = St U Mk [t/2) S v U Mk (t/2) . (6.16) 

Originating in the non-commutativity of Hj^M k an d HjU k , f° r k ^ k' this echo pulse sequence has to be applied for 
each A4k interaction separately and cannot be implemented globally on the total pulse. 

Changing between M.\ and M.2' The addition of the manifold interaction M.\ and M.2 requires to change 
between the configurations Qi = Q Mi and Qi — —Qj Mi G neighb(j). In practice this means that for changing from Q 
to — Q the corresponding phase of the microwave coupling has to be changed by A</> = tt. Alternatively one can also 
keep the laser configuration fixed and apply additional pulses on the quantum states, i.e. 



qA/" C TT Q AT C t TT 



(6.17) 



with S^f° = exp ^— Z7r/2 ^ ieJ ^ c a z z ^j and J\f c the space of non-neighbouring qubits. 



6.1. Creating a cr z cr z -operation by purely global operations? 



In the attempt for only using global interactions on the NV center array, one might think of the conceptually simpler 
method of adding up only Hamiltonians with global unique decoupling fields, i.e. the Hamiltonians B.M\ and Hm 3 

H M k = J2Y at + l £ (6 - 18) 

i i>j 

iGneighb(j) 
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with Mi characterized by Qi = Q Mi and Ms by ^ = — Q Mi. This allows to create a pure a z a z type coupling noting 
that 



H Ml + H Ms — ^2 Jij a i G ] — ^ H zz 

i>j 
iEneighb(j) 

and using e.g. the Suzuki- Trotter pulse sequence (with Um*. = exp(— zHm^)) 



(6.19) 
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(6.20) 



This concept seems to have several advantages: First it requires only global decoupling microwave fields and second it 
does not require to work in the interaction picture (beside the original rotating frame with the laser frequency) thus 
avoiding all the additional complications discussed in the previous section. However there exists a crucial difference 
in the error scaling which now depends on the the magnitude of the decoupling field instead of the much weaker 
dipolar coupling. Therefore timesteps such that Qt/n < 1 are required and thus t/n ~ I/O, ~ r, recalling that an 
effective decoupling is only obtained when the Rabi frequency of the decoupling field is of the order of the inverse 
noise correlation time (see figure S7 a). Such high oscillation frequencies between the decoupling configurations M\ 
and Ms however destroy the decoupling effect as shown in figure S7b with effectively the decoupling field averages 
to zero for fast flipping frequencies. To conclude, creating a o z o z type coupling by the time addition method and at 
the same time preserving the decoupling effect is only possible in the effective interaction picture Hamiltonian frame 
and therefore cannot be performed by purely global addressing. 
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FIG. S7. (a) Fidelity for creating a a z a z -tt/2 pulse for a linear chain of three (blue) and six qubits (red) by directly adding Hmi 
and Hm 3 according to (6.20) with a decoupling field strength Q = 1.2 MHz and a dipolar coupling J = 26 kHz. (b) Oscillating 
decoupling field in the filter spectrum description for t = 25 /us. The red curve illustrates the noise spectrum (T2 = 13.3 /xs, 
r = 2.5 /is, yj (b 2 ) = 30.2 kHz) and the dashed black lines the filter functions for zero decoupling and Qo — 0.5 MHz, 
respectively. The blue, pink and orange lines correspond to filter functions for an oscillating decoupling Q(£) = f2 cos(cj osc £) 
with uJqsc — 1kHz, uJosc — 300 kHz and cj osc = 1MHz, respectively, showing that the decoupling effect vanishes for high 
oscillation frequencies comparable to the noise correlation time. 



6.2. Heisenberg Chain simulation 

As a second application beside the cluster state creation discussed in the main text, we illustrate the possibility 
of simulating a Heisenberg-chain Hamiltonian of the form H — —1/2 E^Li n J v G i with \i = {x,y,z} and 
Jy — Jz — J \ Jx — 8 J- This task essentially requires two steps: Creating the individual parts by using the methods 
described in the previous discussion, namely the <j j x <jj +1 -contribution and the trivial a^ y +1 + a 3 z a^ 1 = S^^ 1 
contributions and in a subsequent step adding those two contribution in time. The accuracy of such a procedure is 
illustrated in figure S8 for a = tt/2 evolution that is limited by non-nearest neighbour interaction as well as the 
number of Trotter cycles (2-3) such that each time slice dt still provides Qdt > 2tt to guarantee the desired H^ k 
Hamiltonian form. 
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FIG. S8. Heisenberg XXZ-chain with 5 qubits and 5 — 1.5. The Hilbert-Schmidt norm distance for a time evolution of 
t — 7r/(2J) is plotted vs the decoupling field and one (blue), two (red) and three (green) Suzuki- Trotter cycles. Dashed lines 
correspond to the pure a\ a 3 z contribution. (Hilbert-Schmidt norm distance: tr[([/ ; — Uy(U' — U)] with U the perfect and U' 
the imperfect realisation.) 



7. SYSTEMATIC ERROR COMPENSATION 
7.1. Compensation cycle for the two qubit gate interaction 

The two qubit gate interaction corresponds essentially to a single qubit rotation in a more complicated two qubit 
manifold Mi and .M2, respectively Introducing the (unknown) systematic error e this results, assuming Vt ^> J in 
the following Hamiltonians 



H IMl ~ ^(1 + e) (sf s 2 + s 1 s%) for ft x = ^ 2 



(7.1) 



what can also be written as 



H IMk ^ J -{l + c)cr^ (7.2) 

with cr^ fe the a x operation denned in the manifolds {| + — ), | — +)} and {| ++), | )} for k=l and k=2, respectively 

(see figure 4 (a) in the main text). 

References [37, 38] provide a method to remove the systematic error contribution based on noting that multiples 
of 2 7r pulses with respect to the ideal error-less gate end up with pure e contributions which can, together with the 
property 

c>4> + cr-0 = 2 cos <p cf x (7.3) 

wherein = cos <p a x + sin <\> a y , be used to create a compensation cycle by means of a Suzuki- Trotter time addition of 
the two contributions in (7.3). Restricting this method to the lowest order of the Suzuki- Trotter expansion (in order 
to keep the pulse sequence simple and additionally because higher orders require more time and therefore end up in 
a loss of additional fidelity by decoherence), it follows that for 2tt cos</> = — 0/2 

Ml e} (7r) M^ } (2tt) M^iir) M* c >(0) = M< e= °}(#) + 0(e 3 ) (7.4) 

with 

M [ ; } {e) = ex V (-i 6 -{l + e)o < /\ (7.5) 
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and 



COS <7 r + sin a*. 



(7.6) 



Instead of giving a detailed derivation of the compensation sequence here we refer to reference [37] and will discuss 
the similar idea for the multiparticle case in 7.2 as well as the special case of the M^it) M^^^M^n) sequence 
in section 7.3 with an emphasis on the complication with increasing particle number. Note that the change of the 
rotation axis <j) can be accomplished by using 



e -i9/2a ct) = e - 



12 a z e -i9/2a x e %4>l2a z 



(7.7) 



what allows to rewrite the rotation as 



M { ;\e) = T^Ml%{6)Tl 



(7.8) 



with = exp(— i(j>/2 a z ). 

The application of (7.4) to the evolution of Hamiltonian (7.2) is straightforward by just replacing 
appearing in (7.5-7.8) by the manifold counterparts <j^ fe , <j^ fc , <j^ fc , <j^ fc . More precise it follows that 



5 &x •> &y i & ' z 



1/2 (ol-ol) 

_1 



^M 2 



1/2 {<tI + <tD 

1 



(J, 



(7.9) 



wherein the last two options have additional contributions that however do not affect the states in the gate manifolds. 
In summary, a compensated gate is obtained using the sequence (see figure S9) 



M| e} (tt) M$ (2 tt) M] e} (tt) M< £ > (<9) = Ml £=0 > (6) + 0(e 3 ) 



(7.10) 



with 



M&(P) = exp(-iH IMh 0/J) 



(7.11) 



2tt cos^ = -0/2 (7.12) 
and M'(0) = T M< c >(0)Tj with T = exp(-z<£/2 a^). 

7.2. Compensated a z a z interaction for multiple qubits 

Assuming systematic errors in the dipole dipole coupling term resulting in the Hamiltonian 

X£ = \ E (l + *«K'i (7-13) 

i>j 
iEneighb(j) 

leads, beside the desired contribution, to additional dependent terms in the defective evolution operator 

M&{0) =ex P X>^j ex P ^ E £ M^°ij • ( ? ' 14 ) 

Noting that the individual contributions in the sum of the evolution operator commute such that it can be also be 
written as a product of the individual components, it is advantageous to analyze first simpler situation 

M iA°) = ex P ( ~ l 7i (1 + e ij) °l a i) = ex P M \ °\ a l J ex P (~i \ e ij °z °l\ • (7-15) 
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FIG. S9. 7r/2 gate in the Mi manifold including systematic error compensation sequence for two qubits. The systematic error 
is characterized by e = 0.2 and the dipolar coupling by J = 26 kHz. No decoherence effects are taken into account. 



A contribution only involving the systematic error part can be obtained by any multiple of 2 7r-pulses 

M£(n27r) = (-l) n exp (-i nn eij <j\ g{) (7.16) 

with n G Z. This property forms the basis for creating a compensation sequence provided that the relatively fixed 
phase (limited by the 2tt condition) can be controlled. The latter is achieved using the properties 

(4 + a%) ai = T</, (4 4) T\ + T_ (<r* a{) = 2 cos cj>4 a{ (7.17) 

with = exp(— i <j)/2 <j l x ) and = coscj)<j z — sine/) a y analogue to (7.6). Performing a Suzuki- Trotter addition of the 
two components cr^ and <r_0 therefore leads to (using (7.16) and (7.17)) 

Mf M (n2n)M? jt _t(n4n)M^^^ (7.18) 

wherein ^(9) — Mfj(0)T^. Comparing this result to the systematic error contribution in (7.15) it is obvious 
that the error part can be compensated by (7.18) if 

Q 

4 n 7r cos = — - (7-19) 
in which case a corrected gate operation is obtained using the sequence 




(7.20) 



Taking into account that the error scales as G(ne^) the parameter n should be as small as possible, leading to the 
obvious choice n = 1. It is straightforward to extend the compensation analysis to the complete evolution (7.14) what 
finally leads to 

Mj e} (2 k)M { _ € 1 (4 tt) Mj e} (2 tt) M< e > (0) = M^ * (0) + 0[e% ) (7.21) 

with M^{0) = T§ M{6) Tj, 8tt cos = — and = exp ^— i (j)/2 ^2^ ie j^ c erf ^ and J\f c the space of non-neighbouring 

qubits. As an example the method is illustrated for the creation of a four qubit cluster state in figure S 10 (a), showing 
that high decoupling fields are required in presence of decoherence in order to benefit from the compensation despite 
the significant increased gate time that is problematic in terms of decoherence. 
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FIG. S10. Multiqubit error compensation (a) Fidelity vs systematic error e (ei2 = e, £23 = — e, 634 = 0.75e) for creating 
a four qubit cluster state without (blue) and with a compensation sequence (red) for Q = 5 MHz (continuous), Q = 10 MHz 
(dashed-dotted) and Q = 1MHz (inset) assuming only nearest neighbour interactions (see section 7.4). (Noise parameters: 
t — 2.5 • 10" 6 s,b = 30.2 kHz, T2q= = 13.3 /xs. Dipolar Coupling: — J — 26 kHz.) (b) Fidelity before (blue) and after 
(red) systematic error compensation (7.21) for a one dimensional cluster state creation of four qubits and assuming an exact 
a z crj creation (neglecting decoherence) . The green dashed lines show the results when applying the advanced (7.22) sequence, 
clearly illustrating the failing of that sequence in the multiparticle case. 



7.3. Failing of the advanced method for the multiparticle case 



As indicated in (7.20) smaller values for n lead to an improved scaling of the compensation method, i.e. to a 
smaller remaining error contribution and in the case when decoherence is relevant, to a significant reduction of the 
compensated gate time. Thus one could ask if also n = 1/2 would be possible. Following the general setup it is 
obvious that for this choice relation (7.16) is not valid any more and therefore e dependent contributions appear in 
the compensation sequences (7.20) and (7.21). However for the case of a single particle rotation, and so also for the 
compensated two qubit gate discussion in section 7.1, there exists exactly such a compensation sequence given by (7.4), 
(7-6) 

M { ; } (ir)M£(2n)M { ; } (ir) (7.22) 

such that for 2tt cos</> = — 0/2 

Mj e} (7r) M^ } (2tt) Afj c} (7r) M< € >(0) = M< €=o >(0) + G(e 3 ) . (7.23) 

Note that, as expected, the choice of rotation angles <ft differs from the one in (7.20). Before analyzing the failing of 
this advanced scheme for the multiparticle case it is useful to consider first its working mechanism for a single particle. 
In that case M < e >(0) = exp (-i 0/2 a z ) and = T rj| with T = exp(-z 0/2 a x ) and accordingly 

Mj e} (7r)M 3 { ^ } (2 7r)Mj e} (7r) = e- i7r/2e ^ \ e -t*/*°+ e -i*e°*+ e -W2^J e -*,r/2 e ^ ( 7>24 ) 

Note that not only e dependent contributions survive due to the fact that two of the contributions in (7.22) are not 
multiples of 2tt. From (7.17), (7.20) and the Suzuki- Trotter expansion it turns out that a valid compensation sequence 
can be constructed if the term in brackets fulfils the following relation 

e — in/2a c j ) g — i tt e cr 3 ^ — in/2cr ( f ) _ ^ — inea-^ 

such that 

Mj e} (tt) Mj e} (2 tt) Mj e} (tt) = e~ l w/2 ea * e~ l 71 6 a -* e~ l 71/2 

v ■ ^ ) 

= exp (—i 2 tt cos (j)ecr z ) . 
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It therefore remains to prove the validity of (7.25). Using relation (7.7) allows to express it as 



;tt/2c 



-i7r/2(j0 -i2cf)cr x p — int 



i2<f>(T x —iir/2cT < f ) 



-( 



i tt/2c 



e 

i 7r/2 cr_0 



(7.27) 



where again condition (7.7) was used together with the crucial requirement that 



e — %7r/2a ( p e ~i(f)<j x _ e +i(f) a x ^ — 1^/2(7^ 



(7.28) 



what follows directly by noting that exp (—i7r/2cr ( f ) ) = —icr^ and {cr^cr^} = 0. 

Now let us return to the case of multiple particles. Herein the same concepts can be applied as in section 7.2 despite 
now using the pulse sequence (7.22). However it turns out that in such a situation the condition analogue to (7.28) 



T<h Si> i |ieneighb(j) < < Tl 



i.g. 



(7.29) 



with and M c defined as below eqn. (7.21), is not fulfilled. That this is in general not valid can best be seen by 
considering the example of three particles in a linear chain in which case 



« *V2 [*l * 3 Z ] e -i °l = e ~i °l e ~i */2 [a\ cj%+cj% a 3 z ] 



(7.30) 



what reveals that due to the double commutation with two cr^ contributions there is no net sign change in the phase 
when commuting exp (— i (j)(J 2 ) in contrast to the situation for a single particle in (7.28). That prevents the application 
of relation (7.7) and hence terms cannot be cast in as in (7.27). One should note that this issue can still be 
fixed for the three particle case by choosing in (7.29) such that the phase change is performed on the first and third 
particle in what case the left and right hand side of that equation would be equal thus allowing to use a compensation 
cycle of the form (7.23). For more than three particles such a possibility however does not exist. 
To conclude: Due to an even number of permutations in (7.29) (i.g. two permutations for a linear chain, four for 
a two dimensional array and six for a three dimensional one) it is in general not possible to extend the n = 1/2 
compensation sequence (7.23) to the multiparticle case (see figure S10(b)). Therefore the general concept (7.21) has 
to be used in that situation, despite being less efficient in both time and e compensation range. 



7.4. Non-next nearest neighbour couplings in the compensation sequence 

The Hamiltonians (7.1) and (7.13) used to describe the compensation mechanism have been defined up to next- 
nearest neighbours, i.e. higher order contributions were neglected so far. However, whereas those contributions 
give only a small correction for simple tt/2 multiqubit pulses, as can be seen from the numerical results in the 
main text, they are crucial in the compensation sequence. Due to the long pulse times of such a sequence higher 
order contributions cannot be neglected. One has to distinguish between two types of higher order couplings: Even 
couplings defined as next-nearest neighbour and third, fifth, . . . nearest neighbour couplings and odd couplings as 
second, fourth, . . . nearest and diagonal couplings. Whereas the first ones have always the form of a cr z a z coupling 
and are automatically compensated for by the next-nearest neighbour based compensation sequence, the second class 
consists just of the reduced manifold Mk couplings and is not removed by the compensation sequence. This failure 
for odd couplings is as well originating from the fact that the compensation pulse acts on both qubits involved such 
that (7.17) is not fulfilled any more, as on the fact that the manifold interaction does not commute with the g z g z 
type contributions such that a splitting of the evolution in a perfect and defective part as in (7.14) is not possible any 
more. 

A solution to that problem consists of eliminating odd order couplings from the beginning by adding different 
contributions in time. Examples for that are given for a two qubit state in figure 5 (b) in the main text and in 
figure SI 1 for second nearest neighbour couplings. With the evolutions obtained that way as a starting point higher 
orders can be removed in a concatenated way taking benefit of the different evolution timescales dependent on the 
qubit distance. For the multiqubit compensation sequence (7.20) it is sufficient to remove odd order couplings up to 
the second order (e.g. up to couplings of qubit one to five), because higher orders effects are too weak for giving 
significant contributions. 
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FIG. Sll. Removing 2nd nearest neighbour interactions for a linear four qubit configuration. Each of the blocks removes 2nd 
nearest neighbour couplings; however both blocks are needed to restore the proper next nearest neighbour interaction. Blue 
lines denote a z a z and red lines Mi interactions and dashed lines denote a negative sign of the corresponding interaction. Red 
marked qubits denote that the interaction is embedded between a U and pulse with U = exp(— in/2a x ). Couplings higher 
than second order are neglected in the illustration. 
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